The R pacakge enerscape computes the energy landscape based on models predicting the energy costs of transport. Two models are currently implemented, the ARC model for terrestrial, legged animals from (Pontzer 2016) and the model for cycling expenditures from (Prampero et al. 1979). Adding other models is relatively easy and can be done by writing a custom function in the enerscape_internals.R script. Currently, there is no implementation for user-specified models, but I plan to add this in the future.
An extensive derivation of the ARC model can be found in Pontzer (2016). Here, it is sufficient to say that the ARC model accounts for the energy spent for cross-bridge cycling, generating tension, as well as associated activation-relaxation processes that restore cellular ion gradients - mostly Ca\(^{2+}\) necessary for cross-bridge cycling. Notably, predictions of the ARC model fit very well observed data. In enerscape, the dataset from Pontzer (2016) can be access as the variable pontzer:
library(tidyverse)
library(raster)
library(sf)
library(enerscape)
#' Function to plot the areas
#'
#' @param x is a raster
#' @param poly is a shapefile polygon
plot_area <- function(x, poly, mask = NULL, col = NULL, void = FALSE) {
if (is.null(col)) {
plot(x, axes = !void, box = !void)
} else {
plot(x, col = col, axes = !void, box = !void)
}
if (!is.null(mask)) {
plot(mask, add = T, legend = FALSE,
col = adjustcolor("white", alpha.f = 0.5),
axes = !void,
box = !void)
}
plot(poly$geometry, add = TRUE)
}
#' ARC model predictions
#'
#' @param mass body mass of the animal (kg)
#' @param slope incline of the terrain (degree)
ARC <- function(mass, slope) {
E_ar <- 8 * mass^(-0.34)
E_mec <- 50 * (1 + sin((2 * slope - 74) / 180 * pi)) * mass^(-0.12)
E <- E_ar + E_mec
return(E)
}
pontzer %>%
as_tibble %>%
mutate(ARC = ARC(Mass, Incline)) %>%
ggplot() +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
geom_point(aes(ARC, Cost.of.Transport, col = Incline)) +
scale_color_gradient(low = "steelblue", high = "tomato") +
scale_x_log10(n.breaks = 11) +
scale_y_log10(n.breaks = 11) +
xlab("ARC predicitons") +
ylab(expression("E"["COT"]*" (J kg"^"-1"*"m"^"-1"*")")) +
theme_bw()To compute energy landscapes using enerscape, we need an elevation rasterLayer, normally called a digital elevation model (DEM). The raster should have equal resolution on both horizontal and vertical axes and a planar coordinate reference systems in meters. The following is a raster that have these characteristics:
## [1] 100 100
## CRS arguments:
## +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs
We can then compute the energy landscape of the area for a general animal of 50 kg using the enerscape() function:
## - Raster cells are assumed to have same horizontal and vertical resolution and with planar coordinate reference system (e.g. UTM)
## | Calculating slope
## | Calculating work
## | Calculating conductance (1 / work)
## - Do not use slope with negative values for distance calculations
The output is an enerscape object, which is a list with elements:
The function en_extrapolation() can be used to highlight where ARC predictions extrapolate from the testing data. A rasterLayer masking extrapolations for the incline is returned if available:
## $`Mass extrapolated`
## [1] FALSE
##
## $`Slope extrapolated`
## [1] TRUE
##
## $`Slope extrapolation`
## class : RasterLayer
## dimensions : 367, 509, 186803 (nrow, ncol, ncell)
## resolution : 100, 100 (x, y)
## extent : 852140, 903040, 4659980, 4696680 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs
## source : memory
## names : Slope
## values : 1, 1 (min, max)
plot(dem)
plot(extr$`Slope extrapolation`,
add = TRUE,
col = adjustcolor("blue", alpha.f = 0.3),
legend = FALSE)
arrows(888523, 4674275, 878791, 4671549)The algorithm to compute least-cost paths from the R package gdistance is implemented in the en_lcp() function (Etten 2017). en_lcp() takes as input two points for which the least-cost path need to be evaluated.1 Alternatively, random points are generated if simulate_random_points = TRUE:
## 1 of 10 random points ...
## 2 of 10 random points ...
## 3 of 10 random points ...
## 4 of 10 random points ...
## 5 of 10 random points ...
## 6 of 10 random points ...
## 7 of 10 random points ...
## 8 of 10 random points ...
## 9 of 10 random points ...
## 10 of 10 random points ...
The random walk algorithm from gdistance is also implemented as en_passage():
## 1 of 10random points ...
## 2 of 10random points ...
## 3 of 10random points ...
## 4 of 10random points ...
## 5 of 10random points ...
## 6 of 10random points ...
## 7 of 10random points ...
## 8 of 10random points ...
## 9 of 10random points ...
## 10 of 10random points ...
## $Origins
## x y
## [1,] 881090 4694830
## [2,] 883990 4676830
## [3,] 870990 4679430
## [4,] 857490 4690330
## [5,] 897890 4688530
## [6,] 876390 4692730
## [7,] 880990 4660130
## [8,] 882590 4692530
## [9,] 886190 4671730
## [10,] 868990 4679330
##
## $Destinations
## x y
## [1,] 854590 4672030
## [2,] 874290 4679530
## [3,] 874990 4692030
## [4,] 880490 4677330
## [5,] 883590 4669330
## [6,] 867490 4661530
## [7,] 893390 4676730
## [8,] 868090 4689830
## [9,] 893390 4671230
## [10,] 869490 4662730
##
## $Passages
## class : RasterStack
## dimensions : 367, 509, 186803, 10 (nrow, ncol, ncell, nlayers)
## resolution : 100, 100 (x, y)
## extent : 852140, 903040, 4659980, 4696680 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs
## names : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7, layer.8, layer.9, layer.10
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## max values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
##
##
## $`Cumulative net passage`
## class : RasterLayer
## dimensions : 367, 509, 186803 (nrow, ncol, ncell)
## resolution : 100, 100 (x, y)
## extent : 852140, 903040, 4659980, 4696680 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer
## values : NA, NA (min, max)
As this may be quite slow, I implemented another way of computing the overall connectivity between two points - or for the whole landscape - using Circuitscape and Omniscape. The initialization files for these to algorithms can be written using circuitscape_skeleton() and omniscape_skeleton(), which create circuitscape.ini and omniscape.ini to be launched from within a Julia terminal. See the example [Landscape connectivity for the Marsican bear (Ursus arctos marsicanus) in the Sirente-Velino Regional Park: Circuitscape and Omniscape implementations.] for details.
Here, I assess the overall landscape connectivity for the Mariscan brown bear (Ursus arctos marsicanus) in the Sirente-Velino Regional Park (SVRP) in central Italy parcosirentevelino.it. In particular, I derive the energy landscape for the area and use it as a “resistance matrix” for animal movement, i.e. assuming that the animal would move across low energy landscape trajectories in order to minimize energy costs.
The required input for enerscape using the ARC model are:
I load the DEM and the park polygon data:
## Reading layer `SVRP-polygon' from data source `/home/eb97ziwi/Proj/enerscape-paper/data/SVRP-polygon.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 30 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 857140.5 ymin: 4665004 xmax: 897975.7 ymax: 4691681
## CRS: unknown
m <- mask(dem, park, inverse = TRUE)
plot_area(dem, park, m, void = TRUE)
title("Elevation (m)")
scalebar(10000,
type = "bar",
divs = 2,
label = as.character(c(0, 5, 10)),
below = "kilometers")I calculate the energy landscape for a bear of 140 kg:
## - Raster cells are assumed to have same horizontal and vertical resolution and with planar coordinate reference system (e.g. UTM)
## | Calculating slope
## | Calculating work
## | Calculating conductance (1 / work)
## - Do not use slope with negative values for distance calculations
Or plotting only the cost of travel per cell:
plot_area(en$rasters$Work, park, m, col = topo.colors(100), void = TRUE)
title("Energy cost of travel (kcal)")I choose two points for which to compute the least-cost path and the connectivity using Circuitscape:
p <- data.frame(x = c(877367, 882653),
y = c(4674192, 4677413))
plot_area(dem, park, m, void = TRUE)
points(p, pch = 20)And I run Circuitscape in Julia.
julia> using Omniscape
julia> run_omniscape("/home/eb97ziwi/Proj/enerscape-paper/data/circuitscape.ini")## [ Info: 2020-12-19 08:31:24 : Precision used: Double
## [ Info: 2020-12-19 08:31:27 : Reading maps
## [ Info: 2020-12-19 08:31:28 : Resistance/Conductance map has 186803 nodes
## [ Info: 2020-12-19 08:31:36 : Solver used: AMG accelerated by CG
## [ Info: 2020-12-19 08:31:36 : Graph has 186803 nodes, 2 focal points and 1 connected components
## [ Info: 2020-12-19 08:31:36 : Total number of pair solves = 1
## [ Info: 2020-12-19 08:31:37 : Time taken to construct preconditioner = 1.683233535 seconds
## [ Info: 2020-12-19 08:31:37 : Time taken to construct local nodemap = 0.008903604 seconds
## [ Info: 2020-12-19 08:31:37 : Solving pair 1 of 1
## [ Info: 2020-12-19 08:31:38 : Time taken to solve linear system = 0.526888286 seconds
## [ Info: 2020-12-19 08:31:39 : Time taken to write current maps = 0.559248897 seconds
## [ Info: 2020-12-19 08:31:39 : Time taken to complete job = 15.111858446
## 3×3 Array{Float64,2}:
## 0.0 1.0 2.0
## 1.0 0.0 31.8625
## 2.0 31.8625 0.0
Sometimes, few cells with extremely high current/connectivity can make the output map hard to read, as most of the connectivity values are squeezed into a narrow range. To avoid this, the value of extremely high outliers can be forced to an upper limit. For example, I force outliers to the 0.999 quantile value of the connectivity and standardize the Circuitscape output between 0 and 1.
cs <- raster("../data/circuitscape_cum_curmap.asc")
cs <- log10(cs)
cs <- cs - min(values(cs), na.rm = TRUE)
cs[cs >= quantile(cs, 0.999)] <- quantile(cs, 0.999) #remove far outliers
cs <- cs / max(values(cs), na.rm = TRUE)
plot_area(cs,
park,
m,
col = viridis::inferno(100),
void = TRUE)
points(p, pch = 3)
lines(lcp$Paths)
title("Connectivity between locations")To compute the overall landscape connectivity, I use the Omniscape algorithm, which runs Circuitscape iteratively on a moving window. I initialize the omniscape.ini file to be passed to Julia:
omniscape_skeleton(en,
file.path("/home",
"eb97ziwi",
"Proj",
"enerscape-paper",
"data"),
radius = 10)And I run the Omniscape algorithm in a Julia terminal. Don’t forget to specify the number of parallel threads: https://docs.circuitscape.org/Omniscape.jl/stable/usage/#Parallel-Processing.
julia> using Omniscape
julia> run_omniscape("/home/eb97ziwi/Proj/enerscape-paper/data/omniscape.ini")## Starting up Omniscape. Using 7 workers in parallel. Using double precision...
## Solving moving window targets...
## Progress: 100% | Time: 0:03:14
## Done!
## Time taken to complete job: 196.3408 seconds
## Outputs written to
## /home/eb97ziwi//home/eb97ziwi/Proj/enerscape-paper/data/omniscape_2
## (Union{Missing, Float64}[0.5152914765395912 0.5884433266779532 … 0.719651357069539 0.8600629588479131;
## 0.4834003129933446 0.6594116484239914 … 0.635221929581483 0.7110678186022065;
## … ;
## 0.2127925084620419 0.2067469779201847 … 0.39577453680538544 0.46681040575465904;
## 0.17382917047741175 0.2093726458198999 … 0.31763823508216504 0.3738227777731409]
## Union{Missing, Float64}[1.0885252650071 1.136985249629082 … 1.169180462521136 1.1104937306261533;
## 1.053217843455199 1.0090282247541849 … 0.9301240063409788 1.1859951520336032;
## … ;
## 1.1798282079383584 0.9871029439774996 … 0.9873244151514602 1.2432973368603717;
## 1.1120398150491893 1.168088759875787 … 1.0583987239714665 1.0877633126993709])
I remove outliers and standardize the Omniscape output between 0 and 1.
os <- raster("../data/omniscape_2/cum_currmap.tif")
os <- os - min(values(os), na.rm = TRUE)
os[os >= quantile(os, 0.999)] <- quantile(os, 0.999) #remove far outliers
os <- os / max(values(os), na.rm = TRUE)
plot_area(os,
park,
m,
col = viridis::inferno(100),
void = TRUE)
title("Landscape connectivity")Here, I study how body weight (kg) and cycling speed (km/h) influence the energy required by participants to cycle L’Eroica event. In particular, I derive the energy landscape for the route and sum the energy costs of all cells.
I load the DEM and the route polygon data:
dem <- raster("../data/eroica/eroica.tif")
route <- st_read("../data/eroica/route.shp") %>%
st_transform(crs(dem))## Reading layer `route' from data source `/home/eb97ziwi/Proj/enerscape-paper/data/eroica/route.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 12 fields
## geometry type: LINESTRING
## dimension: XYZ
## bbox: xmin: 11.3464 ymin: 43.03679 xmax: 11.61873 ymax: 43.49132
## z_range: zmin: 129.4 zmax: 620
## geographic CRS: WGS 84
r <- st_coordinates(route)[, 1:2] %>%
Line() %>%
Lines("Eroica")
r <- SpatialLines(list(r))
crs(r) <- crs(route)
plot_area(dem, route, void = TRUE)
scalebar(d = 10000,
type = "bar",
below = "kilometers",
label = as.character(c(0, 5000, 10000) / 1000),
div = 4)I force the DEM on the route only:
The energy landscape for a cyclist of 60 kg cycling at an overall speed of 30 km/h can then be computed using to use the model from (Prampero et al. 1979) by specifying method = "cycling":
## - Raster cells are assumed to have same horizontal and vertical resolution and with planar coordinate reference system (e.g. UTM)
## | Calculating slope
## | Calculating work
## | Calculating conductance (1 / work)
## - Do not use slope with negative values for distance calculations
The total energy costs of cycling L’Eroica are then derived by summing the energy expenditures for the whole route:
## [1] 5782.049
To get the energy costs for cyclists of different weights and for other cycling speed, I make a for loop, e.g.:
library(foreach)
res <- foreach (weight = seq(50, 100, by = 10),
.combine = "rbind") %do% {
res <- foreach (speed = seq(10, 50, by = 5),
.combine = "rbind") %dopar% {
en <- enerscape(el,
weight,
"kcal",
4,
method = "cycling",
v = speed)
cbind(speed,
sum(values(mask(en$rasters$Work, el)),
na.rm = TRUE))
}
cbind(weight, res)
} %>%
as_tibble()I can then plot the energy costs for different body weight and cycling speed:
res %>%
mutate(weight = as.factor(weight)) %>%
ggplot() +
geom_line(aes(speed, V3, col = weight)) +
geom_violin(aes(speed, V3, group = as.factor(speed)), alpha = 0, draw_quantiles = 0.5) +
geom_point(aes(speed, V3, col = weight), size = 1) +
xlab("Speed (km / h)") +
ylab("Energy consumed (kcal)") +
scale_color_brewer("Weight (kg)", palette = "RdGy") +
theme_classic() +
theme(panel.grid.major.y = element_line(color = "grey80", size = 0.25))and the difference in minimum and maximum costs at same speed, but different weight:
res %>%
group_by(speed) %>%
summarize(Min = min(V3),
Max = max(V3)) %>%
mutate(Diff = Max - Min) %>%
ggplot() +
aes(speed, Diff) +
geom_smooth(alpha = 0.2, col = "tomato") +
xlab("Speed (km / h)") +
ylab("Maximum energy costs\ndifference (kcal)") +
geom_point() +
theme_classic()## `summarise()` ungrouping output (override with `.groups` argument)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Etten, Jacob van. 2017. “R Package Gdistance: Distances and Routes on Geographical Grids.” Journal of Statistical Software 76 (1): 1–21. https://doi.org/10.18637/jss.v076.i13.
Pontzer, Herman. 2016. “A Unified Theory for the Energy Cost of Legged Locomotion.” Biology Letters 12 (2): 20150935. https://doi.org/10.1098/rsbl.2015.0935.
Prampero, P. E. di, G. Cortili, P. Mognoni, and F. Saibene. 1979. “Equation of Motion of a Cyclist.” Journal of Applied Physiology 47 (1): 201–6. https://doi.org/10.1152/jappl.1979.47.1.201.
A path can also be drawn on map using en_path() specifying draw = TRUE and the number of points to draw the path, e.g. (en, draw = TRUE, n = 10).↩︎